Waterlike thermodynamic anomalies in a repulsive-step potential system 
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We report a computer-simulation study of the equilibrium phase diagram of a three-dimensional 
system of particles with a repulsive step potential. The phase diagram is obtained using free-energy 
calculations. At low temperatures, we observe a number of distinct crystal phases. We show that 
at certain values of the potential parameters the system exhibits the water-like thermodynamic 
anomalies: density anomaly and diffusion anomaly. The anomalies disappear with increasing the 
repulsive step width: their locations move to the region inside the crystalline phase. 

PACS numbers: 61.20.Gy, 61.20.Ne, 64.60.Kw 



Some liquids (for example, water, silica, silicon, car- 
bon, and phosphorus) show anomalous behavior in the 
vicinity of their freezing lines 0, [3, & EL [3, @, ■ The 
water phase diagrams have regions where a thermal ex- 
pansion coefficient is negative (density anomaly), a self- 
diffusivity increases upon pressuring (diffusion anomaly), 
and the structural order of the system decreases upon 
compression (structural anomaly) [3, 0. The regions 
where these anomalies take place form nested domains in 
the density-temperature [3] (or pressure-temperature 0) 
planes: the density anomaly region is inside the diffusion 
anomaly domain, and both of these anomalous regions 
are inside the broader structurally anomalous region. In 
the case of water these anomalies are usually related to 
the anisotropy of the intermolecular potential. However, 
isotropic potentials are also able to produce density and 
diffusion anomalies, ft is interesting that such potentials 
may be purely repulsive and can be considered as the sim- 
plest models for the water-type anomalies. It has been 
shown that water-like structural, thermodynamic, and 
dynamic anomalies can be generated in systems where 
particles interact via isotropic potentials with two char- 
acteristic length scales, with shorter range correspond- 
ing to a hard-corelike steep repulsion and longer range 
representing softer repulsion - potentials in which two 
preferable interparticle distances compete depen ding on 
the thermodynamic conditions of t he sy stem pl ,[ 3 ,JlCl l .llll. 

13, [13, [13, EHl, E3, E3, EE [S3, EH, Ha [13, [3 ■ in 

these studies was found that there is an exception case 
- the repulsive-step potential - in which no anomalies 
were reported yet [27[. In this sense, it is very interest- 
ing to mention the recent study [28| of evolution of the 
behavior of the water-like anomalies in the system with a 



tunable potential ranging from a ramp potential, which 
has all mentioned above anomalies to the repulsive-step 
potential, where no anomalies were found so far [23]. In 
[281 ] it was shown that potentials in which two preferred 
distances are present always exhibit water-like anoma- 
lies, but sometimes they are in an inaccessible region, as 
inside a crystal phase. This is the case for the repulsive- 
step potential studied in Ref. [23] ■ 

However, recently it was shown that water-like anoma- 
lies can exist in the systems of particles interacting 
through the repulsive step potential [13] for some values 
of the potential parameters. 

This potential was introduced in the early work of 
Hemmer and Stell [3, [3] in order to describe isostructural 
phase transitions in materials such as Ce or Cs and is the 
simplest example of a repulsive intermolecular potential 
that has a region of negative curvature in the repulsive 
part, a feature that is known to be present in the inter- 
atomic potentials of some pure metallic systems, metal- 
lic mixtures, electrolytes and colloidal systems. Systems 
of particles interacting through such pair potentials can 
possess a rich variety of phase transitions and thermo- 
dynamic anomalies, including liquid-liquid phase transi- 
tions 
region 



30 L 3^ 32 L and isostructural transitions in the solid 

~|iaHl|. 



In this sense, the purpose of this paper is straightfor- 
ward. We will show that the water-like anomalies do 
exist for the repulsive step potential, but with increasing 
the width of the repulsive step they move to the inacces- 
sible region inside the crystal phase. The width_of the 
repulsive step of the potential considered in Refs. 
corresponds exactly to this limiting case. 
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The repulsive step potential has the form: 



$(r) 




(1) 



where d is the diameter of the hard core, a is the width 
of the repulsive step, and e its height. In the low- 
temperature limit T = ksT/e << 1 the system reduces 
to a hard-sphere systems with hard-sphere diameter a, 
whilst in the limit T >> 1 the system reduces to a hard- 
sphere model with a hard-sphere diameter d. For this rea- 
son, melting at high and low temperatures follows simply 
from the hard-sphere melting curve P — cT/a' 3 , where 
c « 12 and a is the relevant hard-sphere diameter (a and 
d, respectively). A changeover from the low-T to high- 
T melting behavior should occur for T — 0(1). The 
precise form of the phase diagram depends on the ratio 
s = a/d. For large enough values of s one should expect 
to observe in the resulting melting curve a maximum that 
should disappear as s — > 1 35]. The phase behavior in 
the crossover region may be very complex, as shown in 

In our simulations we have used a smoothed version 
of the repulsive step potential (Eq. ([1])), which has the 
form: 



$(r) 



" ' +-e(l-tanh(fc (r-<r s ))) (2) 



where n — 14, ko = 10. We have considered the follow- 
ing values of a s : a s = 1.15, 1.35, 1.55, 1.8. In Fig. [I] the 
repulsive step potential is shown along with its smooth 
version which was used in our Monte-Carlo (MC) and 
molecular dynamics (MD) simulations. 

In the remainder of this paper we use the dimensionless 
quantities: f = r/d, P = Pd 3 /e, V = V/Nd 3 = 1/p. As 
we will only use these reduced variables, we omit the 
tildes. 

In [2^] the phase diagrams of the repulsive step poten- 
tial system were reported for a s = 1.15, 1.35, 1.55. In the 
present article we also calculate the phase diagram of the 
system for a s — 1.8. To determine the phase diagram at 
non-zero temperature, we performed constant-NVT MD 
simulations combined with free-energy calculations. In 
all cases, periodic boundary conditions were used. The 
number of particles varied between 250, 500 and 864. No 
system-size dependence of the results was observed. The 
system was equilibrated for 5 x 10 6 MD time steps. Data 
were subsequently collected during 3 x 10 6 5t where the 
time step St = 5 x 10~ 5 . 

In order to map out the phase diagram of the system, 
we computed its Helmholtz free energy using the thermo- 
dynamic integration: the free energy of the liquid phase 
was computed via thermodynamic integration from the 
dilute gas limit [36| , and the free energy of the solid phase 
was computed by thermodynamic integration to an Ein- 
stein crystal [36j, [37 1. In the MC simulations of solid 
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FIG. 1: A repulsive step potential consisting of a hard core 
plus a finite shoulder (dashed line) (e = 1, a = 1.5) along 
with the continuous version of the potential ([2]) used in the 
simulations (e = l,a s — 1.55). 



phases, data were collected during 5 x 10 4 cycles after 
equilibration. To improve the statistics (and to check for 
internal consistency) the free energy of the solid was com- 
puted at many dozens of different state-points and fitted 
to multinomial function. The fitting function we used 
is a p . q T p V q , where T and V — 1/p are the temperature 
and specific volume and powers p and q are connected 
through p + q = N. The value N we used for the most 
of calculations is 5. For the low-density FCC phase N 
was taken equal to 4, since we had less data points. The 
transition points were determined by a double-tangent 
construction. 

The region where we have expected thermodynamic 
anomalies is situated close to the glassy phase, that 
means that proper sampling of the phase space can be 
problematic. To overcome this problem we have used 
the parallel tempering method [36j . Instead of simulat- 
ing one system we consider n systems, each running in 
the NVT ensemble at a different temperature. Systems 
at high temperatures go easily over potential barriers and 
systems at low temperatures sample the local free energy 
minima. The idea of parallel tempering is to put over 
MD the MC scheme of accepting/rejecting a move, but 
in our case it would be accepting or rejecting a swap of 
temperatures between different configurations after each 
full (equilibration together with sampling) MD run. If 
the low and high temperatures are far apart, the prob- 
ability to exchange the configurations is quite low, that 
is why we use a range of 'intermediate' temperatures be- 
tween them with a small temperature step. So after run- 
ning the whole parallel tempering scheme we get a row 
of systems with subsequent temperatures and each of the 
systems was sampled several times. For our problem we 
usually used 8 temperatures and tried to swap them 40 
times. This simulation took almost 24 hours running it 
on 8 processors in parallel at the Joint Supercomputing 
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Center of Russian Academy of Sciences. 

Fig. [2] shows the phase diagrams that we obtain from 
the free-energy calculations for four different values of a s 
(we included the phase diagrams for a s = 1.15, 1.35, 1.55 
for completeness). Fig. HJa) shows the phase diagram 
of the system with a s — 1.15. One can see that for 
the system with a s = 1.15 there are no maxima in the 
melting curve. In a soft-sphere system described by the 
potential 1/r 14 a face-centered cubic crystal structure has 
been reported 38]. However, the addition of a small 
repulsive step leads to the appearance of the FCC-BCC 
transition shown in Figs. Eta). 

Fig. Eljb) shows the phase diagram of the system with 
cr s = 1.35 in the p — T plane. There is a clear maxi- 
mum in the melting curve at low densities. The phase 
diagram consists in two isostructural FCC parts corre- 
sponding to close packing of the small and large spheres 
separated by a sequence of structural phase transitions. 
This phase diagram was discussed in detail in our pre- 
vious publication [29( ]. It is important to mention that 
there is a region of the phase diagram where we have 
not found any stable crystal phase. We think that no 
crystal structure is stable in this de nsity range because 
of frustration as it was discussed in [29J|. In [29j it was 
shown that the glass transition occurs in this region with 
T g = 0.079 at p = 0.53. The apparent glass-transition 
temperature is above the melting point of the low-density 
FCC and FCO phases (see Fig.Etb)). This suggests that 
the "glassy" phase that we observe is thermodynamically 
stable. This is rather unusual for one-component liq- 
uids. In simulations, glassy behavior is usually observed 
in metastable mixtures, where crystal nucleation is ki- 
netically suppressed. One could argue that, in the glassy 
region, the present system behaves like a "quasi-binary" 
mixture of spheres with diameters d and <j s and that the 
freezing-point depression is analogous to that expected in 
a binary system with a eutectic point: there are some val- 
ues of the diameter ratio such that crystalline structures 
are strongly unfavorable and the glassy phase is stable 
even for very low temperatures. The glassy behavior in 
the reentrant liquid disappears at higher temperatures. 

One can expect the frustration to be even more pro- 
nounced if we increase the step size. In Fig. Etc) we show 
the phase diagram of the system with the potential (E| 
for a s — 1.55. One can see that the system also demon- 
strates low- and high density FCC phases separated by 
FCC to BCC transitions and the amorphous gap which is 
much more wider than for a s = 1.35. We did not find any 
crystal structure between these isostuctural phases in our 
study. The glass transition temperature is T g = 0.11091 
at p = 0.5. One can see that the glassification tempera- 
ture becomes higher. Given the lack of crystal structure 
between crystalline phases and the increase of the glass 
transition temperature one can assume that the frustra- 
tion effects become higher with the increase of the step 
width. 
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FIG. 2: Phase diagram of the system of particles interact- 
ing through the potential (2) with <j s — 1.15, 1.35, 1.55, 1.8 
in p — T plane. In Figs. 2 (b-d) it is shown the behavior of 
the diffusivity as a function density. In Figs. 2 (b-c) we also 
represent the locations of the minima on the isochores (see 

Fig.rj) 
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The phase diagram of the system with a s = 1.8 is 
shown in Fig.[2jd). One can see that inside the disordered 
gap in the phase diagram there appears the crystalline 
phase with diamond structure, however, this phase does 
not extend over the whole disordered region in the phase 
diagram. 

As it was mentioned above, one can expect the appear- 
ance of thermodynamic anomalies in the vicinity of the 
anomalous points on the phase diagrams of the repulsive- 
step potential system. To check this point, we calculated 
the isochores and diffusivity for different values of a s . 
In this sense, it is not surprising that there are no ther- 
modynamic anomalies for a s — 1.15. It is known, that 



for normal liquids the diffusivity decreases monotonically 
with increasing density at constant temperature. In con- 
trast, we have observed in our model, that for a certain 
values of the potential parameters, for the densities in the 
vicinity and above the maximum of the melting curve, the 
diffusivity curve has a bend (see Figs. E^b-c)). In Figs. [3] 
the behavior of the diffusivity is shown in more detail for 
different values of a s . One can see that with increasing 
the width of the repulsive step a s the anomaly is becom- 
ing less pronounced and disappears for a s — 1.8. It is 
interesting to note that this value of a s corresponds to 
width of the repulsive step considered in 27|, [28| where 
no anomalies were found for the repulsive step potential. 
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FIG. 5: Re-scaled part of the phase diagram with locations 
of the minima on the isochores. 



The region where the diffusivity anomaly exists almost 
coincides with a region in which the isochore has a min- 
imum instead of growing monotonically (see Figs. [2b- 
c), where the locations of the minima of isochores are 
shown, and Figs. Using the thermodynamic relation 
(dP/dT) v = ap/Kx, where ap is a thermal expansion 
coefficient and is the isothermal compressibility and 
taking into account that Kj- is always positive and finite 
for systems in equilibrium not at a critical point, one 
can conclude that there is a range of densities and tem- 
peratures where the thermal expansion coefficient ap is 
negative. 

To elucidate the behavior of the anomalies with in- 
creasing the width of the repulsive step of the potential, 
we re-scaled the parts of the phase diagrams correspond- 
ing to the first maximum on the melting curve (see Fig. [2j) 
by multiplying the density by the erf and dividing the 
temperature by T max , where T max is the temperature 
corresponding to the maximum. In accordance with the 
qualitative picture depicted above the re-scaled parts of 
the phase diagrams should coincide. As it is seen in Fig. [5] 
this is approximately the case for a s = 1.35,1.55,1.8. 
The discrepancies between the curves appear to be be- 
cause we consider the smoothed version of the repulsive 
step potential. In Fig. [5] we also show the locations of 
the minima of the isochores for a s = 1.35, 1.55. One can 
see that with increasing the width of the repulsive step 
the line of the minima moves to the melting line and be- 
comes invisible in the metastable region. It should be 
noticed that this scenario is similar to the one depicted 
in Ref . S 

At low densities, we have effectively a liquid consisting 
of spheres with diameter er s , at high densities, the liquid 
consists of spheres with diameter d. In the "anomalous 
region" inbetween, our system appears as a mixture of 
both sorts of particles, and one can expect that in this 
region structural order should decrease for intermediate 
values of a s . In this case, the entropy of the system 
should increase with increasing density, and, due to the 



thermodynamic relation (^r) = P 2 (§p) (if) & 
one gets the anomalous behavior in this region. This fur- 
ther demonstrates that our model shows a quasi-binary 
behavior. 

In summary, we have performed the extensive com- 
puter simulations of the phase behavior of systems de- 
scribed by the soft, purely repulsive step potential © in 
three dimensions. We find a surprisingly complex phase 
behavior. We argue that the evolution of the phase di- 
agram may be qualitatively understood by considering 
this one-component system as a quasi-binary mixture of 
large and small spheres. Interestingly, the phase dia- 
gram includes two crystalline FCC domains separated 
by a sequence of the structural phase transitions and a 
reentrant liquid that becomes amorphous at low temper- 
atures. The water-like anomalies (density anomaly and 
diffusion anomaly) were found in the reentrant liquid for 
a s = 1.35, 1.55. The anomalies disappear with increas- 
ing the repulsive step width: their locations move to the 
region inside the crystalline phase in the vicinity of the 
maximum on the melting line. 
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